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Abstract 

In these lectures, we introduce finite temperature QCD on the lattice to non- 
experts of the subject. We first formulate lattice QCD both at zero and finite tem- 
peratures. Then a section is devoted to the topic of improved lattice actions which 
are becoming an essential ingredient of precision studies of QCD on the lattice. We 
then discuss about finite temperature SU(3) gauge theory, i.e. QCD without dynamical 
quarks (quenched QCD). Finally, we report recent status of studies in full QCD taking 
into account the effects of dynamical quarks. 



*' Lectures presented at the 1997 Yukawa International Seminar (YKIS'97) on "Non-Perturbative QCD 
— Structure of the QCD Vacuum — " , Kyoto, Japan, 2-12 Dec. 1997. To be published in the proceedings 
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1. Introduction 



Because of the asymptotic freedom in the UV region, quantum chromodynamics (QCD) 
is strongly coupled in the IR region. Therefore, it is difficult to study the vacuum structure 
of QCD by perturbation theory Two characteristic properties of the QCD vacuum are quark 
confinement and the spontaneous breakdown of the chiral symmetry. On the other hand, 
IR properties are affected by temperature. Actually, studies on the lattice have shown that, 
when the temperature becomes sufficiently high, the low temperature hadronic phase of QCD 
turns into the high temperature quark-gluon-plasma phase, in which quarks are liberated 
and chiral symmetry is restored. The quark-gluon-plasma phase is expected to be realized in 
the early Universe and also possibly in heavy-ion collision experiments. A non-perturbative 
study is required to understand the QCD vacuum at finite temperatures. 

In a conventional formulation of quantum field theories, we first have to regularize diver- 
gences appearing in a perturbative expression of loop corrections, and then perform renor- 
malization, order by order in perturbation theory, in order to obtain finite results by removing 
the divergences. This formulation is deeply based on the perturbation theory and can not 
be applied to a non-perturbative study. Therefore, in order to investigate QCD at finite 
temperatures, we need a definition of QCD that does not resort to perturbation theory. 

This leads us to introduce the lattice as a non-perturbative regularization of the theory: 
When we define field variables on a 4-dimensional hyper-cubic lattice with the lattice spacing 
a, the Fourier transforms of lattice fields are periodic in momentum space, so that we can 
restrict all momenta to the first Brillouin zone —n/a < < n/a. Therefore, a finite lattice 
spacing a naturally provides us with an UV cut-off of 0(1/ a). 

When the lattice volume is also finite, the theory is finite and well-defined. Of course, 
the continuum field theory is obtained in the double limits of infinite lattice volume and zero 
lattice spacing. In these limits, the IR and UV divergences are recovered. Before taking these 
limits, however, we can introduce different calculation techniques. Therefore, a procedure 
for a non-perturbative calculation of field theory is as follows: 

1. Formulate the model at finite lattice spacing a and finite lattice volume V. 

2. Perform non-perturbative calculations. 

3. Take the limits a — > and V — > oo. 

For a non-perturbative calculation in the second step, we can perform numerical simulations 
using super-computers, as well as analytic studies using strong coupling expansions etc.. 
Numerous developments and ideas in the last two decades both in algorithms for numerical 
computations and in the computer technology, have made lattice field theory one of the most 
powerful tools to compute the non-perturbative properties of QCD. 
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In these lectures, we attempt to introduce the basic formulation of lattice QCD and its 
applications to finite temperature physics. In Sec. |2], QCD is formulated on the lattice both at 
zero and finite temperatures. In Sec. ^, recent developments in improving the lattice action 
are discussed. Sec. f| is devoted to the results in finite temperature pure gauge theories. 
Finally, in Sec. [5], the recent status of finite temperature QCD simulations with dynamical 
quarks is discussed. A brief summary is given in Sec. ^|. 

There are a few standard text books on lattice gauge theories The development in the 
field is quite fast. The status of lattice QCD is summarized in reviews in the Proceedings of 
recent international symposium on lattice field theory "Lattice XX" S*. 



§2. Formulation of QCD on the lattice 

2.1. Euclidian field theory 

Lattice field theory is based on a path-integral representation of the field theory. 

Z= /[d0]e i5M , (2-1) 



where S[4>] = J d 3 x.dtL((f), <9 M 0) is the action. In order to apply powerful techniques developed 
in statistical mechanics, we consider the theory in the euclidian space-time by substituting 
t with an imaginary time, t — > —ix^, with real X4. 

Z = |[d0]e- 5EM , (2-2) 

where Se = —iS is the euclidian action. In the euclidian space-time, the propagator of a free 
particle is just an exponential function in the coordinate space, where the mass appears as 
inverse correlation length. In the following sections, we drop the suffix E from the euclidian 
action. 

Lattice discretization of the Euclidian space-time in ( p-2| ) defines the lattice field theory 
we shall study. For definiteness, we consider 4-dimensional hyper-cubic lattices, unless stated 
otherwise. The lattice points are called "sites" and the bonds connecting the nearest neighbor 
sites are called "links". 

Matter fields are introduced on the sites. Then, the simplest lattice action can be ob- 
tained by replacing the derivatives in the continuum euclidian action by lattice differenti- 
ations: d^((){x) — ► 7T- [4>{x + ft) — (p(x — ft)] , where fi is the lattice unit vector in the /i-th 
direction with the length a. In the limit a — > 0, the lattice action smoothly recovers the 
continuum action. The objective of a lattice field theory is to compute the continuum limit 
of expectation values: 

{0) = \ I [#] m 6 " 5W ' Z = I [#] 6 ~ m - (2 ' 3) 



3 



2.2. Gauge theory 



In the continuum, a gauge field A^(x) is introduced to intermediate the different local 
gauge transformations at infinitesimally neighboring points x and x + dx^. Therefore, 0^<9 M 0, 
for example, becomes invariant when we substitute <9 M by the covariant derivative that con- 
tains Apix). When the two points are apart by a finite distance, we should consider the 
"connection", defined as a path-ordered integration of along a path connecting these 
points: 



U (x, y) = P.O. exp 



ig i dz^A^z) 

path: x— yy 



(24) 



Under a local gauge transformation <p{x) — > V{x) <j>(x), with V(x) an element of the gauge 
group, the connection transforms as 



U{x,y) — ► V(x)U(x,y)V\y), 
so that 0's at two points can form an invariant by 

(f) j (x)U(x,y)(j)(y). 



(2-5) 



(2-6) 



Therefore, on a lattice where neighboring points are always separated by a finite dis- 
tance, the fundamental variable for gauge degrees of freedom is the connection. The basic 
formulation of lattice gauge theories was invented by Wilson in 1974§. For sufficiently small 
lattice spacing, the connection between the neighboring sites x and x + p, is given by 



U{x,x + p) = U X>IX = exp [igaA^x)] . 



(2-7) 



We consider U Xyfl to reside on the link connecting x and x + jl, and call it the "link variable". 
The link variables have an orientation. A link in opposite direction (the connection from 
x + ji to x) is given by [t4,^- 

In the pure gauge theory, where we have no matter fields acting as sources or sinks, the 
link variables must form closed loops in order to be gauge invariant. This quantity is called 
the "Wilson loop" .0 The simplest loop is a loop along an elementary square with four links, 
which we call the "plaquette". Therefore, the simplest gauge action consists of plaquettes 
only. 



S, 



gauge 



where 



X,flU 



TJ TJ . rrt rrt 



(2- 



(2-9) 



*) In the SU(jV c ) gauge theory, N c link variables can also join at a site by the totally antisymmetric 
tensor e. Closed loops with such joint points are also gauge invariant. 
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is the plaquette at x in the (/i, v) plane. Using the identification (|2-7p and the Hausdorff 
formula e x e y = e x+y+ ' x,y ^ 2+ ", we can show that ( |2-8| ) recovers the continuum gauge action 
/ d Xj^Tr-^F^ in the limit of vanishing a, provided that the coefficient (3 is given by 

(3 = 2N c /g 2 . (2-10) 

Because the link variables are elements of the compact gauge group SU(iV c ), invariant 
integration (Haar integration) over them is finite. Therefore, on a lattice, we are not required 
to introduce gauge fixing. 

2.3. Continuum limit 

By appropriate redefinitions of the fields and the variables by means of the lattice spacing 
a, the lattice field theories can be defined only by dimensionless quantities. Accordingly, on 
a computer, we have no dimensionful quantities at the beginning. 

The scale a is obtained only through a physical interpretation of results for dimensionless 
quantities. A conventional way to fix the scale in QCD is to identify the rho meson correlation 
length £ with the inverse rho meson mass in the lattice units l/m p a. Then the lattice scale 
is given by 

' MeV" 1 . (2-11) 



770 x £ 

Any dimensionful quantity can be used to fix the scale; m p , the charmonium 1S-1P hyperfine- 
splitting, etc.. A conventional choice for the case of pure gauge QCD is the string tension in 
the static quark potential. 

From the relation fl2Tl|) , we see that a smaller a corresponds to a larger £. In particular, 
the limit a — > is achieved at £ — ► oo. We are interested in the modes with correlation 
length of 0(£) in this limit. 

In the coupling parameter space of a lattice theory, £ —>■ oo is achieved at a second 
order phase transition point. Physics of IR modes (modes with large £) near a 2nd order 
phase transition point can be described by the theory of critical phenomena. In statistical 
systems, the critical phenomena have the property of universality, i.e., they do not depend 
on the details of the microscopic theory. We, therefore, expect that the continuum limit of 
a lattice theory is independent on the details of the form of the lattice action. One of many 
consequences of this universality is the recovery of rotational symmetry in the continuum 
limit, because the physics becomes insensitive to the choice of the lattice orientation. 

In QCD, £ — > oo is realized at (3 — > oo (g — > 0) due to the asymptotic freedom. Because 
this is the limit of weak coupling, we can compute the scale dependence by perturbation 
theory. We can formulate the perturbation theory on the lattice similar to the cases in the 
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continuum theory, using the identification (2-7). The result for the scale dependence from 



lattice perturbation theory is given by the beta-function: 

a^- = b g 3 + b l9 5 + --- (2-12) 
da 

where the first two coefficients b = (l/16vr 2 )(ll-2A^ F /3) and b x = (1/16tt 2 ) 2 (102-38iV F /3) 
are universal, i.e. equal to those for the beta-function in the continuum QCD. 

Therefore, in asymptotically free theories, the gauge coupling parameter dependence of 
physical quantities near the continuum limit is under control. This feature is quite important 
to extract precise predictions from the results obtained on lattices with finite a. 

2.4. Fermions 

For fermions on the lattice, a complication exists in the formulation. When we naively 
discretize the Dirac action in the continuum, we obtain the "naive" lattice fermion action: 

S naivc = a 4 E {E ^ I ****'/"* + ™o , (2-13) 

where \P X is a 4-component Grassmann field residing on the site x. Then the propagator 
shows the pole for p 4 = — iE at 



sinh 2 Ea = m 2 a 2 + ^ sin 2 (p*a) (2- 14) 

Near the origin of the momentum space, (|2-14j) predicts the energy eigenvalues to be E 



±ymg + Y^iPi as expected. However, in addition to this mode, we also have 7 extra low 
energy modes inside the first Brillouin zone, at momenta p ~ (ir/a, 0, 0), (0, ir/a, 0), • • •, i.e. 
ir/a for more than one spatial components of p. In the continuum limit, these additional 
modes also contribute to the partition function and the expectation values. In other words, 
they survive as relevant dynamical freedoms in the continuum limit. These unwanted modes 
are called "doublers". 

The problem of doublers is essentially caused by the fact that the derivative in the 
kinetic term of the Dirac action is first order. This causes the term sin 2 (j9ja) in ( 2- IT ) 



to vanish not only at the origin but also also at p — (vr/a, 0, 0) etc. Even if we introduce a 
more complicated lattice derivative than the naive action (|2-13|) , the corresponding lattice 
derivative in the momentum space changes sign from negative to positive near the origin of 
the momentum space. Because the same happens at the origins of all other Brillouin zones, 
as far as we assume analytic continuity, the lattice derivative necessarily crosses zero again 
within the first Brillouin zone. This leads to the additional low energy modes. 
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There exists a rigorous statement, the No-Go theorem by Nielsen and NinomiyaB, that as 
far as we require hermiticity, locality and chiral symmetry, doublers are inevitable. Therefore, 
in order to formulate a lattice theory for one flavor of Dirac fermions in the continuum 
limit, we have to abandon some of these nice properties. The violated properties should be 
recovered in the continuum limit. In the continuum limit, we naively expect that different 
fermion formalisms lead to a universal continuum limit, recovering all the violated properties. 
Off the continuum limit, it is important to check the lattice artifacts due to the formulation 
of lattice fermions. 

There have been a few proposals to formulate fermions on the lattice. Two of them are 
commonly used in major simulations; the Wilson fermion formalism!! and the staggered 
(Kogut-Susskind) fermion formalism!!* . In the Wilson fermion formalism, the chiral symme- 
try is violated on a finite lattice, and in the staggered fermion formalism, locality is violated 
when we try to define a single flavor fermion.0 

2.4.1. Wilson fermion 

In the formulation of the Wilson fermion we introduce a second- derivative term, the 
'Wilson term", to the action: 

4 \- lf# &x+ik + ^x-a - 2W X _ 4 ^ - 1 - cos(p M a) 
a- a }^^x tt^t = a 2^% —% (2-15) 

that is of order a for the mode at the origin p ~ 0, but acts as an 0(1 /a) mass term to 
the doublers. Therefore, the doublers are decoupled in the limit a — > 0, leaving the physical 
mode at the origin untouched. Introducing a dimensionless field ip x = a 3 ^ 2 \P x /\/2K, with 
K = 1/(8 + 2moa), the Wilson fermion action is customarily written as 



^Wilson = ^2i ! x D x<y ijjy (246) 

D x , y = 5 x , y - K^2{(1 - 7 M ) S x+fl) y + (1 + 7m) Ky+p} ( 2 ' 17 ) 



I'- 

The parameter K is called the hopping parameter, and corresponds to the freedom of the 
bare fermion mass: 

1/1 1 \ , 

with K c = 1/8 the point where the bare mass mo vanishes. 
In a gauge theory, the kernel D XyV is modified as 

D x ,y = 5x, y - K {(1 - 7m) U x ,, <W,, + (1 + 7^) Ui 3fl <Wa} > (2-19) 



*) Recently, several new ideas to construct a chiral fermion have been studied. See reviews on the issue 
for details©. 
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where U x ^ is the link variable. 

For Np flavors of quarks, we simply sum up quarks with different flavors. Clearly, the 
flavor symmetry is manifest. However, because the Wilson term is essentially the mass term 
for doublers, the chiral symmetry is violated even in the limit of vanishing bare quark mass. 
As a result, the global flavor/chiral symmetry of continuum QCD, SU(iVir)L X SU(./Vf)r x 
U(l)v, is explicitly broken down to SU(iVV)v x U(l)v- 

Through a perturbative study of axial Ward identities, it is shown that the effects of chiral 
violation due to the Wilson term are removed by appropriate renormalizations, including an 
additive renormalization of the quark mass, i.e. by a shift of K c as a function of 

Away from the perturbative region, we can define K c as the points where the pion 
mass vanishes at zero temperature, or alternatively where the quark mass vanishes at zero 
temperature. Here, the quark mass is defined through an axial-vector Ward identity @>'§, 

2m g Z P {0\P\ir(p= 0) ) = -m n Z A ( | A 4 | n(p = 0) ) (2-20) 

where P is the pseudoscalar density and A± the fourth component of the local axial vector 
current, with Zp and Z& renormalization factors.0 Both definitions give the same K c when 
the chiral symmetry is spontaneously broken 0). K c forms a smooth and monotonic curve 
connecting 1/8 in the weak coupling limit (J3 = oo) and 1/4 in the strong coupling limit 
(f3 = 0). K c can be considered as the points where the chiral symmetry is effectively 
recovered. 

Aoki proposed an alternative interpretation of the massless pion at K c without resorting 
to the chiral symmetry in the continuum limit 0. In this picture, the massless pions are 
identified with the Goldstone modes associated with a second order phase transition, and a 
rich phase structure was predicted at K > K c . Although the physical region relevant to the 
continuum limit of QCD is below the i^-line, understanding the system at K > K c is useful 
in studying the phase structure of Wilson quarks, in particular, at finite temperatures 0-*. 



See also discussions in Refs. 10) and 13) 



2.4.2. Staggered (Kogut-Susskind) fermion 

In the formulation of the staggered fermion, in order to reduce the number of doublers, 
the Grassmann field x on the sites have only one component B>. 

'S'stag = / ] Xx Qx,y Xy (2'21) 



*> We can alternatively define the quark mass by the perturbative formula m q = (Z m /2a) (K^ 1 — A"" 1 ). 
Although these different definitions of m q give different values at finite (3, it is shown that they converge to 
the same value in the continuum limit. EP 
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Qx, y = m a6 XiV + { Uxfi6 



ttW}- ( 2 - 22 ) 



The conventional four component Dirac spinor is constructed collecting the x fields dis- 
tributed on a hypercube. Because a hypercube has 2 4 sites, we end up with 4 flavors of 
degenerate Dirac fermions. 

In this formalism, the flavor-chiral symmetry SU(iVV)L x SU(./Vf)r x U(l) of massless 
iV^-flavor QCD is explicitly broken down to U(A^^/4)l x U(]VV/4)r due to a flavor mixing 
interaction at a > 0. However, at least a part of the chiral symmetry is preserved. Therefore, 
the location of the massless point m = is protected against quantum corrections by this 
symmetry. This makes a numerical analysis of chiral properties much easier than in the case 
of the Wilson fermion for which K c must be determined numerically. It is also known that 
the staggered fermion requires less computer resources (memory etc.) than those for the 
Wilson fermion. 

The action (2-21) can describe quarks only when Np is a multiple of 4. A usual trick 
for the physically interesting cases Np = 2 and 3, is to modify by hand the power of the 
fermionic determinant in the numerical path-integration. 



This necessarily makes the action non-local, which sometimes poses conceptually and tech- 
nically difficult problems. 

2.5. Finite temperature 

Finite temperature field theory is defined by the Matsubara formalism for finite tempera- 
ture statistical systems: We consider static problems in thermal equilibrium at temperature 
T. Instead of the time coordinate, we introduce a coordinate ranging from zero to l/T, 
which formally looks like an euclidian time. (We set ks = 1 in the following.) Then, the 
canonical partition function Z can be written as 



with H the Hamiltonian and S[(f>] = J dx^ J d 3 x L((j), d^cj)). This expression of Z is formally 
equivalent to the path-integral representation of the partition function for an euclidian field 
theory. The differences are (i) the range of the euclidian time 24 is [0, l/T]. and (ii), 




(2-23) 




(2-24) 



according to the trace operation in ( [2-24 ), bosonic and fermionic fields obey periodic and 



anti-periodic boundary conditions in the euclidian time direction. 
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Fig. 1. Speed of dedicated computers developed by lattice physicists, together with those of recent 
commercial machines. 



A lattice discretization of ( |2-24|) defines a finite temperature lattice field theory. In order 
to approximately realize the thermodynamic limit, the lattice size in the spatial direction 
must be sufficiently larger than the size 1/T in the time direction. 

Denoting the dimensionless lattice size in the time direction as N t , the temperature is 
given by 

T = 1/N t a. (2-25) 

In QCD with small Np, a is a decreasing function of j3. Therefore, when Nt is fixed as in 
many numerical simulations, larger {3 corresponds to higher temperature. 

2.6. Numerical simulations 

When we try to perform a numerical evaluation of the euclidian path integral (|2-3|) , we 
quickly encounter the following two difficulties: First, the dimension of integration is huge. 
Even on a quite small lattice, for example, 10 4 , a 10 000 dimensional integration is required for 
each degrees of freedom. Second, the integrand drastically changes its magnitude, forming a 
quite sharp peak in configuration space. This is caused by the factor e~ s in ( j2-3|) . The peaks 
correspond to the classical solutions and the width of the peaks the quantum fluctuations. 
Therefore, a naive mesh method to evaluate an integral cannot give a reliable value until 
the mesh becomes extremely fine. These difficulties can be solved by introducing the Monte 
Carlo method: In order to evaluate the integral, we generate the sample points with a 
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probability proportional to e~ s , in place of the mesh points in the naive method. In the 
text books, we can find several techniques to generate sample points with the correct weight; 
the heat-bath method, the Metropolis method, the Langevin method, etc. The number of 
sample points required for a given accuracy is much smaller than that required in the mesh 
method, and has no strong dependence on the dimension of the integration. 

Nevertheless, when we want to obtain a precise result, which is indispensable in phe- 
nomenological studies, the required computer power is enormous. In order to perform a 
reliable extrapolation of physical quantities to the continuum limit by suppressing lattice 
artifacts, the lattice should be fine enough while keeping a sufficiently big physical volume. 
Simultaneously, we need a high statistical accuracy which requires a large number of sam- 
pling points. 

In order to generate a configuration (a sample point) in the SU(3) pure gauge theory, the 
conventional pseudo heat-bath method© requires about 5 700 floating point operations per 
link. Vector computers in the '80s and parallel computers in the '90s have supported the de- 
velopment of numerical simulations in lattice QCD. As one of the largest user groups in high 
performance computing, lattice physicists have also contributed a lot to the development of 
computer technology itself. Several groups of lattice physicists have even developed parallel 
computers dedicated to lattice QCD©'©©. Fig. 1 shows the speed of recent dedicated 
machines for lattice QCD. On the CP-PACS constructed at the University of Tsukuba0, 
more than 10 000 configurations on a 64 3 x 112 lattice can be generated in a day. 

When the system includes dynamical fermions, different algorithms had to be developed 
as we cannot have Grassmann variables directly on the computers. However, even with the 
latest algorithm several hundred times more computer time is required for fermions. 
Therefore, major QCD simulations on large lattices have been performed in the approxi- 
mation that dynamical pair creation/annihilation of quarks are neglected (quenched QCD). 
Realistic simulations of QCD with dynamical quarks (full QCD) have just begun on recent 
high performance computers, by combining big computer power with the idea of improved 
lattice actions which shall be discussed in Sec. BL 



§3. Improved lattice actions 

3.1. From BETTER to MUST 

In the scaling region near the continuum limit, we expect that the lattice results repro- 
duce the continuum results at distances larger than about the correlation length. However, 
in general, expectation values at short distances in lattice units deviate from the contin- 
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a) b) c) 

Fig. 2. Loops of length 6; (a) rectangle, (b) chair, and (c) parallelogram. 

uum results even at large f3. Therefore, in order to extract a prediction for the continuum 
limit, precise data at large distances on a correspondingly large lattice are required. Such 
simulations are quite expensive, in particular, near the continuum limit. 

On the other hand, due to the universality, there are infinitely many candidates for the 
lattice action describing the same continuum limit, i.e. we can introduce additional terms to 
the standard action without affecting the continuum limit. For example, the lattice gauge 
action may contain non-minimal loops besides the plaquettes adopted in the standard action: 
Using loops up to length 6, we can write 

Sgauge 

-P { c o X p P la q + Cl X 

-Prect + C2 X ^chair + C3 X Pparaj , (3T) 

where P rec t, -Pchair, and P para are rectangle, chair, and parallelogram loops shown in Fig. ||. 
The coefficients in ( |3-1|) should satisfy a relation c + 8ci + 16c 2 + 8c 3 = 1 in order to obtain 
the conventional continuum gauge action in the limit a — > with the identification Q2-10 ). 
Hence, three parameters are left free. From the universality, the long distance behavior of the 
theory is insensitive to these additional parameters. However, short distance properties do 
depend on these parameters. Therefore, a judicious choice of the additional parameters may 
suppress the short distance lattice artifacts even at moderate values of the lattice spacing. 
Such actions are called "improved actions". 

Recently, much progress has been reported in improving lattice actions©'©. When an 
appropriate improved action is obtained, we will be able to perform a reliable extrapolation 
to the continuum limit using data obtained on coarse lattices. It is, of course, better to apply 
such an action to save money on computer time. However, through recent developments in 
numerical studies of lattice QCD, we now have stronger motivations to improve the action. 

As the first motivation, we would like to introduce the recent results of high statistic 
simulation in quenched QCD by the CP-PACS Collaboration!*. The quenched light hadron 
spectrum is computed for Wilson quarks on lattices with the spatial size ~ 3fm to an accu- 
racy of 1-3% in the continuum limit. The baryon spectrum turns out to show a systematic 
disagreement which is of the order 10% (maximally 10 standard deviations) from the exper- 
imental values, which is considered as evidence of systematic errors due to the quenching 
approximation. Therefore, in order to test QCD to an accuracy better than 0(10%), we 
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have to perform full QCD simulations without the quenching approximation. Full QCD sim- 
ulations are, however, extremely computer time consuming compared to those of quenched 
QCD. Even with the TFLOPS-class computers that are becoming available, high statistics 
studies, indispensable for reliable results, will be difficult for lattice sizes exceeding 32 3 x 64. 
Since a physical lattice size of L m 2.5-3.0fm is needed to avoid finite-size effects, the small- 
est lattice spacing one can reasonably reach will be a ~ O.lfm. Hence lattice discretization 
errors have to be controlled with simulations carried out at lattice spacings larger than this 
value. This will be a difficult task with the standard plaquette and Wilson quark actions 
since discretization errors are of order 20-30% even at a « O.lfm. 

We also encounter difficulties with the standard action in a study of finite temperature 
QCD. In the simulation of a finite temperature system, the spatial lattice size should be 
sufficiently large in order to approximately realize the thermodynamic limit. With the 
limitation of the computer power, this means that we have to suppress the lattice size in the 
time direction N t . Currently N t « 4-6 is used in most simulations with dynamical quarks. 
Because the temperature is given by T = 1/N t a, the corresponding lattices are coarse; for 
T ~ T c « 100-200MeV, a ~ 0.2-0.4fm. Subsequent lattice artifacts sometimes make the 
analysis and interpretation of lattice results not a straightforward task. 

The basic idea behind improvement may be obtained by considering the lattice derivative. 
The naive derivative A x f(x) = ^-[f(x+a) —f(x—a)] converges to the continuum derivative 
with errors of 0{a 2 ). We can reduce this error down to 0(a 4 ) by replacing A x — > A x — ^-A x 
which operates on fields at up to next nearest neighbor sites. Similar substitution is effective 
also to obtain a lattice action that approaches to the continuum action much smoothly. Here, 
however, what we want to obtain is not a smoother lattice action, but an action which leads 



to smaller discretization errors in the physical observables (|2-3| ) containing all the quantum 
corrections. 

Several different strategies have been proposed to obtain such a lattice action. Two 
major approaches are the renormalization group (RG) improvement programs first proposed 
by K.G. Wilsonii, and the perturbative improvement programs started by K. Symanzik©. 
In the following subsections, I describe these methods in more detail. Irrespective to the 
differences in approach, the final goal of improvement is to obtain a lattice action which 
shows scaling from a coarser lattice. 

In general, improved lattice actions contain additional terms (interactions that usually 
have a wider spatial size) compared with the standard action. Because such additional 
terms make the simulation quickly difficult and computer time consuming, the efficiency of 
improvement should be tested for each of the new terms introduced. 
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3.2. Symanzik improvement 

In Symanzik's method, we improve physical quantities using a result of perturbative 
computations It consists of the following procedures: 

(i) Compute a set of physical quantities in perturbation theory using a lattice action with 
non-minimal interactions. 

(m) Adjust the non- minimal coupling parameters to remove the leading finite a deviations 
from the continuum limit in these physical quantities. 

A conventional choice for the physical quantities in this program are the low- dimensional 
operators in the effective action. When discretization errors are eliminated up to the nth 
order in a, the resulting action is said to be u O(a n ) improved". For example, the standard 
plaquette gauge action has 0(a 2 ) errors. Several different actions removing these errors, 
to tree or one-loop level, have been proposed depending on the choice of additional terms 

This improvement program is quite attractive because an improved action can be ob- 
tained by an analytic calculation. Nevertheless, the method in its naive form remained 
unsuccessful for a long time. The reason is that the perturbation theory, using the bare 
gauge coupling constant as the expansion parameter, has poor convergence, and the per- 
turbative results do not agree with the results from numerical simulations which are mostly 
performed at (3 = 0(1). 

The failure of the bare perturbation theory may be understood as follows. The bare 
perturbation theory is based on an expansion of the link variable (|2-7|) 



U Xtll « 1 + igaA^x) - \ g 2 a 2 A^xf + ■ ■ ■ . (3-2) 



In order that the perturbation theory works well, the higher order terms in Q3-2| ) should 
be much smaller than 1. However, in actual simulations, the link variable deviates much 
from unity; we obtain the plaquette expectation value about 0.4-0.6. Perturbatively, this 
deviation is caused by large contributions of "tadpole" diagrams©: By fixing the gauge 
appropriately, the mean value of the link variable can be expressed as 

u = (U Xtfi ) wl-i g 2 a 2 {A^xf) + • • • . (3-3) 

Because of a quadratic divergence in the tadpole contribution (A^x) 2 ), higher order terms 
[g 2 a 2 (A^(x) 2 )] k are suppressed only by g 2k , instead of the naive factor g 2k a 2k . In the simu- 
lations, g 2 ~ 0(1). 



14 



3.2.1. Mean-field improvement 

The convergence of the perturbation theory can be improved by an appropriate choice of 
the expansion parameter. Recent significant progress in Symanzik improvement is based on a 
combination of the Symanzik improvement program with an idea of "mean-field improvement 
(tadpole improvement)" of the perturbation theory©'©: Consider a modified expansion of 
the link variable around the mean- field value Uq\ 

U X)fl = u exp iguFaA'^x) . (3-4) 

With much smaller quantum fluctuations, we will obtain a much better converging pertur- 
bation theory with A' (x). From a substitution 

'S'gauge — — P ^ -Pplaq ¥ ~~ U ft ~~ 4 -fplaqj (3' 5) 

plaq plaq ^0 

we find <7mf — 9/ u o- A conventional choice for u$ is uq = (P) 1 / 4 . The perturbation theory 
using guF as the expansion parameter is shown to agree with numerical results much more 
precisely even at the large bare coupling g used in numerical simulations 

Therefore, a prescription for a mean-field improved Symanzik improvement is to introduce 
a third step: 

(Hi) Multiply an appropriate power of ito to the perturbatively determined coefficients of 
the Symanzik action. 

It is, however, not trivial whether the non-perturbative quantities we are interested in are 
also sufficiently improved by this method. Non-perturbative tests are required to confirm 
the efficiency. 

3.2.2. Symanzik-improved Wilson quark action 

The Wilson quark action has 0(a) errors. Symanzik-improved Wilson quark action was 
studied by Sheikholeslami and Wohlert©. The 0(a) improved action reads 

i 

Sclovcr = ^Wilson + - Csw K 1p x O ^ v F x ^ u lf) x , (3'6) 

where a^ v = ^[7^, 7^]. Here F X41V is the field strength on lattice. Because we conventionally 
construct F x ^ v in terms of four plaquette-like loops in the \w plane around the site x, the 
action is called the "clover action" . 

The coefficient cg\y, the clover coefficient, is 1 at the tree-level. Its mean-field improve- 
ment can be done by cgw — ► csw/wjj because the field strength gives the factor Uq 4 while 
the hopping parameter is redefined as K — > K/uq^'®. 
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3.2.3. Non-perturbative determination of a Symanzik action 

Recently an approach to combine the Symanzik method with a non-perturbative study 
was proposed: 

(i') Measure a set of physical quantities by a numerical simulation using a perturbatively 
inspired form of the Symanzik action, varying the coupling parameters freely. 

(ii 1 ) Adjust the non-minimal coupling parameters of the action such that several physical 
requirements are satisfied by the numerical data. 

An attempt to determine csw non-perturbatively requiring a PCAC relation to hold is re- 
ported in Ref. |32[). A technical background for this is the development of the Schrodinger 
functional method in which the lattice correction to the PCAC relation is sensitive to the 
value of csw- 

3.3. RG improvement 

In order to see the basic idea of RG improvement © , let us consider, as a typical example, 
a block transformation of scale factor 2. A correlation function GbiockfV) on the blocked 
lattice is related to the corresponding correlation function G rig(?") on the original lattice by 
Ghiockij) ~ G or i g (2r), where the distance r is measured in lattice units. Therefore, when we 
have the continuum behavior at distances r > ro on the original lattice, then the continuum 
behavior is realized at tq/2 on the blocked lattice. This means that a block transformation 
improves the action. 

A block transformation generally induces many effective interactions in the effective ac- 
tion: 

exp {-Sblockf^block]} = J [d^orig] -fr^block, 0orig] exp { -S^rig (0orig) } , (3"7) 

where the kernel -K^biock, 0ori g ] defines the block transformation orig — > 0biock- Therefore, 
repeated applications of a block transformation defines a flow (RG flow) in the infinite 
dimensional coupling parameter space of effective actions. 

For the SU(3) gauge theory, we imagine the RG flows to be as shown in Fig. |^. In this 



figure, ci, C2, etc. denote additional coupling parameters [cf. the action (|3- . The points 
on the (3 axis correspond to the standard one plaquette action at various (3. Because a 
block transformation makes the correlation £ in lattice units smaller, RG flows are oriented 
towards smaller (3. The hyperplane at (3 = oo (£ = oo) is called the critical surface. When 
the starting action is in the scaling region that is close to the critical surface, we expect 
that all trajectory flows first along the critical surface, and then leaves the critical surface to 
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gradually approach to a universal curve, the "renormalized trajectory" (RT), after sufficiently 
many steps, because the long distance physics is universal in the scaling region. The number 
of block transformations required to get sufficiently close to RT is nothing but log 2 of the 
minimum distance to obtain continuum behavior with S^°\ 

We expect that the RT is also a RG flow starting from an infra-red fixed point, P^, on the 
critical surface. Because the RT is directly connected to the action P^ in the continuum limit 
by block transformations, actions on the RT show continuum properties from the shortest 
distances on the lattice. Therefore, the actions on the RT are called "perfect actions" 

If an infinite number of coupling parameters are admitted, a perfect action is a goal of 
improvement. In reality, we can not handle infinitely many couplings. A few approaches 
are proposed to obtain an improved action in a finite dimensional subspace of the coupling 
parameters. 

3.3.1. Fixed point action approach 

Hasenfratz and Niedermayer proposed to use an approximation of the fixed point action 
Poo as an approximate perfect action©. In asymptotically free theories, because the critical 
surface is located in the weak coupling region, a compact integral equation for P^ can be 
written. With a proper Ansatz for the effective action, we can numerically solve the fixed 
point action P^. The most delicate point is the choice of the finite-dimensional Ansatz action 
("parametrization"). It is noted that generally, a large number of coupling parameters are 
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required to achieve a good approximation for P c 



Many developments are reported in the literature; see Refs. |2(|),[34]) for more details. For 



a trial to compute the RT nonperturbatively, see Ref. The fixed point action for full 
QCD is discussed in Refs. [Zg),|3§,[37p. 



3.3.2. Iwasaki's method 



In an alternative approach, proposed in 1983 by Iwasakit^, the practical constraint for 




the number of coupling parameters in numerical simulations is taken more seriously. In this 
approach, instead of trying to obtain an approximate perfect action directly, one attempts 
to accelerate the approach to the RT by taking advantage of the extended parameter space: 

(i) We restrict ourselves to a small dimensional coupling parameter space consisting only 
of interactions which can be easily implemented in numerical simulations. 

(ii) Find a set of coupling parameters that minimizes the distance to the RT after a few 
block transformations. 

Because a block transformation induces many effective couplings, the number of coupling 
parameters for the initial action can be quite small. 

As an illustration, let us consider an action Si mp in the two-dimwnsional coupling param- 
eter space (/3, Ci) shown in Fig. |, and suppose that the additional parameter c\ is adjusted 
such that the RG flow gets sufficiently close to RT after, say, 2 block transformations. This 
means that, when we simulate the system using Si mp , a separation of 2 2 a is enough to get 
the continuum properties. 

Iwasaki applied this program to the SU(3) gauge theory B 1 . Asymptotic freedom again 
helps us to compute such S imp : 



where Cq = 1 — 8ci and C\ = —0.331 (—0.293) to minimize the distance to RT after one (two) 
block transformations. This action is remarkably simple. It is easy to write a vectorized 
and/or parallelized program for this action. Good efficiency in removing several lattice 
artifacts is reported by the Tsukuba group©'©. An application to finite temperature QCD 
with dynamical quarks will be discussed in Sec. 5.1.2. 

§4. SU (3) gauge theory at finite temperature 

As the first step towards finite temperature QCD, in this section, let us study the case 
of the SU(3) pure gauge theory, or, equivalently, QCD in the approximation that dynamical 




(3-8) 
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creation/annihilation of quarks are neglected (quenched QCD). Although quenched QCD is 
not quite realistic, we can obtain high statistics data on large lattices in this case. Therefore, 
it provides us with a good lesson for full QCD studies. 

4.1. Deconfining transition in SU(3) gauge theory 

In a study of the deconfining transition in pure gauge theories, it is useful to consider 
the "Polyakov loop"© 

^=^Trfn^ )t;4 J. (4-1) 

In the path integral, f2$ is just the factor from the current appearing when a static charge 
is located at the spatial point x. Therefore, except for the renormalization corrections, 
{[}) ~ e~ Fq / T , where F q is the free energy of the static charge. We expect that F q = oo 
(i.e. (f2) = 0) when charge is confined, while F q < oo ((£2) > 0) when charge is deconfined. 
Therefore, (i?) is an order parameter for the deconfining transition. 
The global symmetry behind this order parameter is given by 

jj _^ I z U x>ft if x 4 = and /i = 4 ^ 
1 U X)IX otherwise 

where z is an element of the center group Z(iV c ) of the gauge group SU(iV c ). Because z 
commutates with any element of SU(iV c ), the plaquettes as well as more extended closed 
loops are invariant under the transformation (4-2), i.e. the pure gauge actions ( plf ), (|3-1|). 
etc. are invariant under (4-2). On the other hand, because Q$ crosses the hyper plane x 4 = 
only once, Q$ — > z f2g under (4-2). Therefore, a non- vanishing vacuum expectation value of 
i? implies the spontaneous breakdown of the center Z(3) global symmetry at the deconfining 
transition. 

We may study the nature of the deconfining transition using an effective theory of the 
Polyakov loops in three dimensions. Because the order-disorder phase transition in a three 
dimensional Z(3) spin model (the Potts model) is first order, we may expect that the decon- 
fining transition in the SU(3) gauge theory is also first order, when the interaction in the 
effective theory of £2$ is short ranged©. 

Figure 4 is a result of the Polyakov loop histogram in the complex plane, obtained just 
at the deconfining transition temperature in the SU(3) gauge theory©. The peak at the 
center i? ~ is the contribution of the symmetric phase (low temperature confining phase) 
and the three peaks at non-zero \Q\ with arg Q ~ and ±27r/3 are from the Z(3) broken 
phases. Clear separation of these peaks, as well as the lattice volume dependence shown in 
Fig. 5, implies that the transition is of first order. The order of the transition is confirmed 
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Fig. 4. Polyakov loop histogram in the SU(3) gauge theory at the deconfining transition point 
obtained on a 24 2 x 36 x 4 lattice©. 
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Fig. 5. Polyakov loop histogram for \ f2\ at the deconfining transition point obtained on (a) Nt = 4 
and (b) 6 lattices©. The simulation temperature for the 24 3 x 6 lattice is slightly lower than 
the transition temperature. 



by a precise finite size scaling test©'©. Accordingly, the short range nature of the effecti 
interaction between the Polyakov loops is also shown©. 
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4.2. Transition temperature 

Precisely speaking, we have no real transition on a finite lattice. The transition point 
/3 c (N t ) for each fixed N t can be computed by an extrapolation to the infinite spatial volume 
using a finite size scaling. The value of (3 c {N t ) is translated into physical units by fixing the 
scale at this value of (5. 

In pure gauge theories, the scale is conventionally fixed with the string tension in a static 
quark potential a at zero temperature. 

aa 2 = - lim - \n(W{A)) (4-3) 

A^oo A 

where W(A) = YieedAUe is the Wilson loop with area A. We then obtain T c j ' y/a = 
l/[N ty Jaa 2 (P c (N t ))]B Using recent data© for aa 2 , we obtain T c / ' yfo~ « 0.63(1) in the 
continuum limit (N t — > oo) from the standard action©'©. Results from various im- 
proved actions are T c /y/a~ « 0.62-0.66il ) 'i)'i)'©. Adopting a phenomenological value 
a = (427MeV) 2 from a charmoniun spectrum 0*, we obtain T c ~ 270MeV. 

4.3. Thermodynamic quantities 

In a phenomenological study of the quark-gluon plasma in heavy ion collisions and in 
the early Universe, it is important to evaluate thermodynamic quantities such as the energy 
density e and the pressure p, near the transition temperature of the deconfining phase tran- 
sition. These quantities are defined by derivatives of the partition function with respect to 
the temperature T and the physical volume V of the system 

1 d\nZ r^dlnZ , , 

< = -7«ft- P=T ^ (4 ' 4) 

On a lattice with the size N% x Nt, V and T are given by V = (N s a s ) 3 and T = l/(N t at), 
with a s and a t the lattice spacings in spatial and temporal directions. Because N s and Nt 
are discrete parameters, the partial differentiations in (4-4) are performed by varying a s and 
a t independently. 

Anisotropy on the lattice is introduced by different coupling parameters in temporal 
and spatial directions. For an SU(-/V C ) gauge theory, the standard plaquette action on an 
anisotropic lattice is given by 

S = -f3 s /'„„ - A E Pxm> ( l-">) 



*' Ironically, the largest error comes from aa 2 determined at zero temperature, because an extraction 
of aa 2 contains many delicate fittings. The systematic errors from the choice of the fitting ansatz, fitting 
range, etc. are not fully estimated. 
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where P x ^ u is defined by Q2-9Q . Then the energy density and pressure, renormalized at T = 0, 
are given byS- 1 '© 

eat = 3£ 2 ((P.) ~ (P)o) + fr ((Pt) ~ (P>o)} (4-6) 

((P s )-(P)o) 

where (P s (t)) is the space(time)-like plaquette expectation value and (P)o the plaquette 
expectation value at the same coupling parameters on a zero-temperature lattice. In these 
expressions, the variables £ = a s /a t and a s are chosen to control the lattice spacings. 

Therefore, in order to compute e and p from the results of simulations, the values for the 
derivatives of gauge coupling constants with respect to the anisotropic lattice spacings (the 
anisotropy coefficients) are required. 

The calculation of these anisotropy coefficients in the lowest order perturbation theory 
was done by Karschi!''. [Coefficients needed for the case with dynamical quarks are also 
computed perturbatively in Ref. |54j)-1 However, the bare perturbation theory is not reliable 
for the values of (3 where MC simulations are performed. Accordingly, the perturbative 
coefficients are known to lead to a pathological result of negative pressure at strong couplings 
used in the numerical simulations. In the case of SU(3) gauge theory, the transition is of 
first order. At a first order transition point, we have a finite gap for the energy density, 
the latent heat, but expect no gap for the pressure. It is known that the perturbative 
anisotropy coefficients have the problem of non-vanishing pressure gap at the deconfining 
transition point: Ap/T 4 = -0.32(3) and -0.14(2) on 24 2 x 36 x 4 and 36 2 x 48 x 6 lattices®. 
Therefore, we need non-perturbative values for the anisotropy coefficients. 

We are mainly interested in the values of the anisotropy coefficients for the case of 
isotropic lattices (f3 s = f3 t = (3, i.e. £ = 1) because most simulations are performed in this 
case. At £ = 1, we have (a s df3 s / da s )^ =1 = (a s df3 t / da s )^ =1 = ad/3/da, where adf3/da is the 
beta-function. Furthermore, a combination of the remaining two anisotropy coefficients is 
known to be related to the beta-function©: (d(3 s /d^)^ =l + (df3 t /d^)^ =1 = — (1/2) ad/3 /da. 
Therefore, we need to estimate nonperturbative values of the beta-function ad/3 / da and one 
independent combination of the anisotropy coefficients. 

In connection to this, it is useful to consider a combination e — 3p which is defined via a 
uniform scale transformation: 

a 4 (I d (3 *\ a 4 / d d\, „ 

e - 3p= -VJf [rWrn +3V dv) lnZ = -Vjf [ at da- t + as da~ s ) ln Z 
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Fig. 6. Plaquette history in the SU(3) gauge theory at the deconfining transition point obtained 
on (a) 24 2 x 36 x 4 and (b) 36 2 x 48 x 6 lattices 0. 



3T*N}a^((P. + P t )-2(P) l 
da 



(4- 



at £ = 1. Therefore, this combination depends only on the beta-function ad/3/da. Several 
non-perturbative values for the beta-function are available; from a MC renormalization group 
study 0-*, from a study of the transition temperature 1 , or from the /^-dependence of a 
physical quantity, such as the string tension©. 

As an application, we compute the energy gap Ae, identified with A(e — 3p) because 
we expect Ap = 0, at the deconfining transition in the SU(3) gauge theory. In order to 
determine the expectation values in each phase, we first inspect the "flip-flops" between the 
confining and deconfining phases in Monte Carlo time histories (see Fig. 6 for example), and 
separate the runs into the two phases by the flip-flops. A sufficient number of iterations 
around the flip-flops and around the spikes should be removed to avoid contamination from 
transition stages. Fig. 7 shows an example of expectation values of observables in each 
phase as a function of the number of removed iterations. We can obtain stable results when 
a sufficiently large number of iterations from the transition stages are removed. Note 



*' When we try to define an expectation value for a phase in terms of a cut for the spatial average of 
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(e-3p)/T 4 - 36 2 x48x6 
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Fig. 7. (e — 3p)/T 4 average in each of the confining and deconfining phases in the SU(3) gauge 
theory at the deconfining transition point obtained on a 36 2 x 48 x 6 lattice©, as a function 
of the removed iterations around the phase flip-flops shown in Fig. 6(b). e — 3p is defined by 
(4-8) using the perturbative beta-function. 



that such unambiguous separation of two phases is possible only on large lattices where the 
persistence time of each phase is sufficiently long. Using the beta-function calculated from 
recent string tension data©, we find that A(e-3p)/T* = 2.072(43) and 1.578(42) for N t = 4 
and 6, respectively, with the standard action 0. Using Symanzik improved actions with 2x1 
loops, values of 1.57(12) and 1.40(9) are reported for tree-level and tadpole-improved actions 
on a 32 3 x 4 latticed . More data at larger N t are needed to make a reliable continuum 
extrapolation. 

In order to determine e and p separately, we need one more input as discussed above. 
A non-perturbative determination of a combination of the anisotropy coefficients was at- 
tempted in Refs. |57D,[58|),p9|),|60D using a matching of space-like and time-like Wilson loops 
on anisotropic lattices (the matching method) Alternatively, we can evaluate a non- 
perturbative value of pressure directly from the Monte Carlo data by the integral method^); 
assuming homogeneity of the system, expected for the case of large spatial volume, we ob- 
tain the relation p = —f, where / is the free energy density, / = — y InZ. For the case of 
pure gauge theory with the plaquette action, / can be evaluated by integrating the plaque- 
tte (P) in term of (3 on isotropic lattices, since ■§g\nZ = 6NfN t (P). The resulting value 

the Polyakov loop etc., the result strongly depends on the value of the cut. Therefore, we cannot obtain a 
reliable result in this way. 
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Fig. 8. Anisotropy coefficients (Karsch coefficients) in the SU(3) gauge theory. Dot-dashed lines 
are the results of the perturbation theory H). Filled squares are the results of a matching 
method obtained in Ref. [59]). Filled triangles and thin lines are the results of a matching 
method© combined with the beta- function using a recent string tension data©. Dashed lines 
are the results of the integral method©. Open circles are the results of a new method using 
the transition curve in the coupling parameter space for anisotropic lattices 0. 



of the pressure, in turn, provides us with a non-perturbative estimate of a combination of 
the anisotropy coefficients©'©'©'©. Recently, a new method, using a measurement of 
the finite temperature deconfming transition curve in the lattice coupling parameter space 
extended to anisotropic lattices, was proposed 0. 

Recent results for the anisotropy coefficients in the form 




(k 



(4-9) 



a a :fixed, 
-2, 



a 3 ifixed, £=1 

(Karsch coefficients) with (3 S = 2N c g~ 2 ^~ 1 and fit = 2N c g^ A £, are summarized in Fig. 8. 
Applying these results, we can reanalyze the pressure gap. At the deconfming transition 
point for N t = 6, we find Ap/T A = -0.003(17) using the results of Ref. ||) or -0.040(43) 
using the result of Ref. ^0|). We find that the problem of a non- vanishing pressure gap is 
removed with non-perturbative anisotropy coefficients. 



4.4. Interface tension 
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Fig. 9. (a) Plaquette histogram in the SU(3) gauge theory at the deconfining transition point 
obtained on a 24 2 x 36 x 4 lattices®. See Fig. 6(a) for the corresponding time history. The 
two peaks are separately fitted with two Gaussian distributions, (b) Mix configurations at 
a first order transition point. Above: naive configuration with bubbles. Below: dominant 
configuration for two phase coexistence on a finite box with periodic boundary conditions. 

Table I. Interface tension crj/T^ computed using the histogram method for the data of \ f2\. Large 
spatial volume limit is taken, except for the results at Nt = 3 where the spatial lattice volume 
is 12 3 . 

N t standard Symanzik action© fixed point 

action© tree-level MF improved action© 
~~ 2 0.092(4)^ 

3 0.2434(24) 0.0158(11) 0.0307(8) 

4 0.0295(21) 0.0152(26) 0.0152(20) 0.026(5) 
6 0.0218(33) 



Because the deconfining transition is first order for SU(3), we expect that hot and cold 
phases can coexist just at the transition temperature. When the transition is also of first 
order in full QCD, the value of the surface tension for the interface between two phases is 
important in a study of hadronic bubble formation in the cooling process of the quark gluon 
plasma at the early Universe, heavy ion collisions etc.. To gain experience for the study in 
full QCD, the interface tension has been measured in quenched QCD using various actions. 

Recent numerical computations of the interface tension are based on the "histogram 
method"©: As shown in Fig. 6, the Monte Carlo time history of an observable that is 
sensitive to the transition shows flip-flops between two phases at the transition point. Ac- 
cordingly, when we plot a histogram from the history, we obtain two peaks corresponding 
to the two phases, when the lattice size is sufficiently large. As shown in Fig. 9(a) for the 
case of the plaquette on a 24 2 x 36 x 4 lattice, the valley between the two peaks is in general 
higher than that expected from an overlap of two distributions corresponding to the two 
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peaks. This excess is due to the contribution of transition stages (mixed states) around the 
flip-flops discussed in Sec. 4.3. Naively, a mixed state would look like bubbles of one phase 
floating in the sea of another phase. With positive interface tension, however, the dominant 
contribution will be the configurations with minimum interface area. See Fig. 9(b). Then, 
because the interface area is approximately fixed from the lattice geometry, we can calculate 
the probability to have such a mixed state as a function of the interface tension. Comparing 
the result with the measured probability obtained from the histogram, we can compute the 
value of the interface tension©. 

The actual procedure in lattice QCD is much more complicated. We perform a finite 
size scaling analysis taking into account various corrections including those from capillary 
wave excitations on the interface 0*. Recent results from various actions are summarized in 
Table I. A naive continuum extrapolation using the data for N t — 4 and 6 from the standard 
action gives <J//T C 3 = 0.16(4), which is close to the values obtained with Symanzik improved 
actions at N t = 4©, but slightly smaller than the value from a fixed point action©. 



§5. Finite temperature QCD with dynamical quarks 

Finally, let us study finite temperature QCD with dynamical quarks. As discussed in 
Sec. 2.6, a full QCD simulation requires several hundred times more computational time 
compared with a quenched simulation to achieve corresponding accuracy. Therefore, in 
many cases, results for full QCD are not yet quite quantitative. 

In this lecture, we concentrate on the topics of the order of the finite temperature transi- 
tion in full QCD. With dynamical quarks, the center Z(N C ) transformation (4-2) is no longer 
a symmetry of the action. Therefore, although the Polyakov loop is still a good "indicator" 
of the deconfining transition (it sensitively changes its magnitude at the transition point), 
the group Z(iV c ) is no longer a good guide to study the nature of the QCD transition. 

On the other hand, when quarks are light, we expect that the chiral symmetry, which is 
spontaneously broken in low temperatures, is recovered when the temperature is sufficiently 
large.0 Using the universality hypothesis, the nature of the finite temperature QCD transi- 
tion near the chiral limit can be studied by a Ginzburg-Landau effective theory respecting 
the chiral symmetry of QCD, the effective a model. From a study of the effective a model 
at finite temperatures©, the transition in the chiral limit (the chiral transition) is predicted 
to depend quite sensitively on the number of light quark flavors Np. Let us consider QCD 

*' Numerical simulations show that, when the quarks are in the fundamental representation of SU(iV e ), 
the deconfining transition in the heavy quark mass limit smoothly turns into the chiral transition when we 
decrease the quark mass. 
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with Np degenerate light quarks. For Np > 3, a first order transition is predicted from 
the sigma model.0 For Np = 2, on the other hand, the order of the transition is not quite 
definite in the effective a model; a first order transition is predicted when the anomalous 
axial Ua(1) symmetry is effectively restored at the transition temperature, while a second 
order transition is expected otherwise. However, because the Ua(1) breaking operator is a 
relevant operator whose coefficient grows towards the IR limit under a renormalization group 
transformation, the transition is more likely to be second order 0). A non-perturbative study 
is required to determine the order of the transition for Np = 2 conclusively. 

In nature, we know six flavors of quarks; u, d, s, c, b, and t. The lightest two quarks, u 
and d, are much lighter than the relevant energy scale for thermal processes near the critical 
temperature T c ~ 100-200 MeV: m u , <C T c . On the other hand, the last three quarks, c, 
b, and t, are sufficiently heavy that they are expected to play no appreciable roles in thermal 
processes near T c . In the following, the case Np = 2 corresponds to the case where the third 
quark s is much heavier than the relevant energy scale; m s ^> T c , while the case Np > 3 
corresponds to the case m s <C T c . Because m s ~ 150-200 MeV is just of the same order of 
magnitude as the expected values of T c , in order to make a reliable prediction for the real 
world, we have to fine-tune the value of m s in the more realistic case of two light quarks and 
one heavy quark (Np = 2 + 1). 

5.1. Chiral transition for Np = 2 QCD 

Understanding the nature of the QCD transition for Np = 2 is an important step toward 
the clarification of the transition in the real world. When the transition in the chiral limit 
(the chiral transition) is second order, we expect that the transition turns into an analytic 
crossover at non-zero m q , while when the chiral transition is first order, it will remain to be 
first order for small m q . In order to confirm the expected crossover numerically, we have to 
study the lattice size dependence to see if the formation of a singularity (e.g. the increase of 
the peak height of a susceptibility with increasing the lattice volume) stops on sufficiently 
large lattices. However, it is difficult to numerically distinguish between a very weak first 
order transition and a crossover, especially at small m q . 

Here, the universality provides us with useful scaling relations that can be confronted 
with numerical results of QCD, in order to test the nature of the transition: It is plausible 
from an effective a model that, when the chiral transition is of second order, QCD with 
two flavors belongs to the same universality class as the three dimensional 0(4) Heisenberg 
model©. The 0(4) model is much simpler than the a model, and its scaling properties 

*) We restrict ourselves to the case Np < 6. See Ref. fol) for the phase structure at Np > 7. 
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Fig. 10. (a) Peak height of the magnetic susceptibility for Np = 2 QCD with staggered quarks 
as a function of the spatial lattice volume Nf. (b) The same data as a function of the lattice 
volume rescaled by zero-temperature pion correlation length. 



are well studied. For example, at small external field h near the critical temperature T c 
for h — 0, the pseudo-critical temperature T pc (h) and the peak height of the magnetic 
and thermal susceptibilities follow T pc — T c ~ h Zg , ~ h~ Zm , and xT ax ~ h~ zt , where 

z g = 1/(35, z m = 1 — 1/5, and z t = (1 — (3)/ ' [35 in terms of the 0(4) critical exponents [3 
and 5. Here the values of (3 and 5 for the 0(4) model are well established©. In QCD, we 
identify T ~ 6/g 2 , h ~ m q , and M ~ (WW). 

In lattice QCD, an additional complication should be noted because no known lattice 
fermions have the full chiral symmetry on finite lattices, as discussed in Sec. ||. In partic- 
ular, on a coarse lattice used in a finite temperature simulation, we sometimes encounter 
sizable deviations from the scaling behavior expected in the continuum limit. Therefore, the 
appearance of the 0(4) scaling is also a useful touchstone to test the recovery of the chiral 
symmetry on the lattice when the chiral transition is of second order. 

5.1.1. Results with staggered quarks 

The 0(4) scaling was first tested on the lattice for staggered quarks by the Bielefeld 
group©. Based on simulations on an 8 3 x 4 lattice at m q a = 0.02, 0.0375, and 0.075 using 
the standard action, they obtained z g = 0.77(14), z m = 0.79(4), and z t = 0.65(7), where 
the corresponding 0(4) values© are 0.537(7), 0.794(1), and 0.331(7). The result for z m 
is consistent with the 0(4) value while other exponents are in disagreement with the 0(4) 
values. 

Possible causes of the discrepancy are (i) m q is not small enough to see the critical 
behavior in the chiral limit, and (ii) the spatial lattice volume is not large enough to obtain 
the observables in the thermodynamic limit. Two additional caveats are in order for Np = 2 
staggered quarks: (iii) The symmetry in the chiral limit at a > is 0(2) instead of 0(4). 
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Practically, however, the values of the 0(2) exponents are almost indistinguishable from the 
0(4) values with the present numerical accuracy, (iv) The action is not local. Therefore, 
an assumption behind the universality argument can be violated so that some non-universal 
behavior may appear©. The correct continuum chiral limit with the 0(4) symmetry will be 
obtained only when we first take the continuum limit a — > and then take the chiral limit. In 
addition to these points, we also have to check technical details in the numerical simulation; 
the accuracy of the methods to simulate the system, such as the the finite step-size error 
and the dependence on the convergence criterion for fermion matrix inversion. 

A systematic study of the quark mass dependence as well as the lattice volume dependence 
is in order. The JLQCD Collaboration performed a series of simulations on 8 3 x 4, 12 3 x 4, 
and 16 3 x 4 lattices at m q a = 0.01, 0.02, 0.0375, and O.O750 } . The Bielefeld group also 
extended their study to larger spatial lattices©. The results obtained are consistent with 
each other. It turned out that determination of critical exponents on 8 3 x 4 lattices suffers 
from a sizable finite lattice-size effect for m q a < 0.0375. 

From the lattice-size dependence of the magnetic susceptibility Xm for a fixed value of 
quark mass, shown in Fig. 10(a), we find that the transition is a crossover for m q a > 0.02; 
the peak height x™ ax f° r m q a = 0-02 stabilizes on spatial lattices larger than 12 3 . For 
m q a = 0.01, on the other hand, x™ ax is increasing up to the largest spatial lattice of 16 3 . 
If this increase is maintained up to infinite volume, then the transition is first order at this 
quark mass. However, no clear indication of a first-order transition are found from the 
lattice volume dependence of Monte Carlo time histories and histograms at m q a = 0.01. 
Furthermore, when the lattice volume is rescaled by the zero-temperature pion correlation 
length, the lattice volume 16 3 for m q a = 0.01 approximately corresponds to the volume 
12 3 for m q a = 0.02, where the increase of x™ ax terminates [see Fig. 10(b)]. Therefore, it is 
possible that the increase of x™ ax f° r m q a = 0-01 * s a transient effect. 

Assuming that the finite size effect is sufficiently small on the 16 3 x 4 lattice, we fit the 
data at the four values of m q a. We find z g = 0.64(5), z m = 1.03(9), and z t = 0.82(12). 
(Removing the data for m q a = 0.01 gives slightly smaller but consistent values with larger 
errors.) The results for z g and z m deviate sizably from the 0(2) or 0(4) values. On the other 
hand, the identity z g + z m — z t = 1 expected for a second-order fixed point with two relevant 
operators is approximately satisfied. Thus, the exponents are consistent with a second-order 
transition at m q = 0. 

In summary, we find that the determination of the nature of the two-flavor chiral transi- 
tion with staggered quarks using the standard action to involve subtle problems. While the 
data so far do not contradict a second-order transition at m q = 0, the exponents take quite 
unexpected values, at least in the range m q a > 0.01. Evidently further work, possibly on 
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Fig. 11. The Polyakov loop and the pion screening mass obtained with Wilson quarks with a RG 
improved action©. In these figures, the conventional notation (5 = 6/g 2 is used. Larger (3 
corresponds to a higher temperature. The horizontal axis 1/K — 1/Kq is proportional to m q a. 



larger spatial sizes and smaller quark masses, is needed to clarify this important problem. 

5.1.2. Results with Wilson quarks 

Let us now study the issue using Wilson quarks. It turned out that Wilson quarks in the 
standard action lead to several unexpected phenomena on lattices with N t = 4 and 6: On 
these lattices, the transition becomes once very sharp when m q is increased from the chiral 
limit contrary to the expectation in the continuum limit that the chiral transition 

becomes weaker with larger m q . Together with other strange behavior of physical quantities 
near the transition point, this phenomenon is identified as an effect of lattice artifacts^. 

Therefore, the Tsukuba group applied an improved action©. Their action is the RG im- 
proved gauge action ( |3-8| ) with C\ = —0.331, coupled with the standard Wilson quark action. 
Although the quark part is not improved, the lattice artifacts observed with the standard 
action are shown to be well remove They also find that the physical quantities are 

quite smooth around the transition point at m q > 0, as shown in Fig. 11. The straight line 
envelop of m 2 at finite temperature (Nt = 4) shown in Fig. 11(b) agrees with m 2 obtained 
at low temperature (N t = 8), and corresponds to the PCAC relation m 2 oc m q expected in 
the low-temperature phase. The smoothness of the physical observables strongly suggests 
that the transition is a crossover at m q > 0. 

Concerning the nature of the transition in the chiral limit, the transition becomes mono- 
tonically weaker with increasing 6/g 2 (see Fig. 11). Because the transition point shifts to 
larger 6/g 2 at larger m q , increasing 6/g 2 corresponds to increasing m q for the transition. In 
the chiral limit m 2 decreases monotonically to zero as the temperature is decreased from 
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Fig. 12. (a) Chiral condensate as a function of h = 2m q a for Wilson quarks with a RG improved 
actionlliP. (b) Best fit for the scaling function with 0(4) exponents with x 2 /df = 0.61, using 
the identification h = 2m q a and t = (3 — j3 c t, where (3 = 6/g 2 and Pct is the value of 6/g 2 
at the chiral transition point. Do not confuse with the critical exponent (3 appearing in the 
combination t/h 1 ^ 6 . The plot contains all data of (a) within the range < 2m q a < 0.8 and 
6/g 2 < 2.0. Solid curve represents the scaling function obtained in an 0(4) spin model. 

above, towards the chiral transition point. At the transition temperature for finite m q , m 2 
also shows a similar monotonic decrease with decreasing m q . These results suggest that the 
chiral transition is continuous. 

For a more decisive test about the nature of the transition, a scaling study is required. 
From the universality argument, the magnetization M in a spin model can be described by 
a single scaling function near a second order transition point: 

M/h 1/s = fit/h 1 ^ 5 ), (5-1) 

where h is the external magnetic field and t = [T — T c ]/T c the reduced temperature. When 
the QCD transition is of second order in the chiral limit, the chiral condensate should satisfy 
this scaling relation with 0(4) exponents 1/(35 = 0.537(7) and 1/5 = 0.2061(9)0 and also 
with the 0(4) scaling function f(x). 

The magnetization M is identified with the chiral condensate in QCD. For Wilson quarks, 
the naive definition of (ipip) f° r the chiral condensate is not adequate because the chiral 
symmetry is explicitly broken due to the Wilson term. A proper subtraction and a renor- 
malization are required to obtain the correct continuum limit. A properly subtracted {ipip) 
can be defined via an axial Ward i dentityS); 

$V)sub = 2m 9 aZ]T(7r(x)7r(0)> (5-2) 

X 

where Z is the renormalization coefficient. For our purposes, it is enough to use the tree 
value, Z = (2K) 2 . When the chiral symmetry is spontaneously broken, a singulariry in the 
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Fig. 13. (a) Map of expected nature of the QCD transition for Np = 2 + 1 QCD as a function of 
the u and d quark mass m u d and the s quark mass m s . (b) The same figure with the results 
from Wilson quarks using the standard action© . First order signals are observed at the points 
marked by filled circles, while no clear two-state signals are found at the points represented by 
the open circles. The values of quark mass in physical units are computed using a determined 
by m p (T = 0) = 770 MeV. The real world corresponds to the point marked by the star. 

pion propagator cancels the factor m q in the r.h.s. of (|5-2j), giving a finite value for (V^sub 
in the chiral limit. 

The results of M = (ipip) sn b around the deconfming transition/crossover for N t = 4 is 
shown in Fig. 12(a). Fig. 12(b) shows the result from a fit of M to the scaling function 
obtained for an 0(4) model©, by adjusting the chiral transition point (3 ct and the scales for 
t and h, with the exponents fixed to the 0(4) values. The scaling ansatz works remarkably 
well with the 0(4) exponents. A recent study shows that the situation holds also when data 
at t < are included 0*. On the other hand, a change of the exponents quickly makes the fit 
worse: For example, fixing the exponents to the MF values, suggested by Kocic and Kogut 
as a possibility for two-flavor qcdB, the data no longer falls on the MF scaling function. 

The success of this scaling test with the 0(4) exponents suggests strongly that the chiral 
transition is of second order in the continuum limit. It also indicates that the chiral violation 
due to the Wilson fermion action is sufficiently small with this improved action, for the values 
of m q and 6/g 2 studied here. 

5.2. Influence of the strange quark 

In order to study the nature of the transition in the real world, we should study the in- 
fluence of the s quark. Our expectation about the nature of the finite temperature transition 
as a function of quark masses is summarized in Fig. 13(a), neglecting the mass difference 
among u and d quarks (Np = 2 + 1). The limit m s = 00 corresponds to the case Np = 2 
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discussed in Sect. 5.1 where we found second order transition at m U d = 0. For m U d = m s 
(Np = 3), the transition is of first order in the chiral limit. Therefore, on the axis m u d = 0, 
we have a tricritical point m* where the second order transition at large m s turns into first 
order il*. For m s < m*, the second order edge of the first order transition region is sug- 
gested© to deviate from the vertical axis according to m u d oc (m* — m s ) 5 / 2 . Our main goal 
of investigations with the s quark is to determine the position of the physical point in this 
map. 

With staggered quarks, Brown et al. 15 found, for the degenerate Np = 3 case (m u d = 
m s = m q ), a first order signal at m q a = 0.025, (3 = 5.132. For Np = 2 + 1, they obtained a 
time history suggesting a crossover for m^a = 0.025, m s a = 0.1 at /3 = 5.171. Their study 
of hadron spectrum at this simulation point on a zero-temperature lattice leads to a m^/m p 
smaller than the experimental value, suggesting that this m s is smaller than the physical 
value. At the same time, their large m^jm p suggests that their m u d is larger than the 
physical value. This implies that the physical point is located in the crossover region unless 
the second order transition line, which has a sharp m u d dependence near m* (cf. Fig. 13(a)), 
crosses between the physical point and the simulation point. In order to obtain a more 
decisive conclusion, a study to systematically investigate a wider region of the parameter 
space in Fig. 13(a) is required. 

The Tsukuba group studied the issue with Wilson quarks using the standard action 0). 
They found first order signals for m q < 140 MeV, while no clear two state signals were 
observed for m q > 250 MeV, where the physical s quark mass, giving = 1.02 GeV, is about 
150 MeV in their normalization of m q . (The scale was fixed by m p at zero tmeperature.) 
For Np = 2 + 1, first order signals are observed for m u d ~ at both m s ~ 150 and 400 MeV. 
A study of zero-temperature hadron spectroscopy for Np = 2 + 1 shows that ~ 1.03(5) 
GeV at the simulation point m s ~ 150 MeV, verifying that this simulation point is very 
close to the physical point. The results are summarized in Fig. 13(b). The physical point is 
located in the first order region. 

Although both staggered and Wilson simulations give a phase structure qualitatively 
consistent with Fig. 13(a), Wilson quarks tend to give larger values for critical quark masses 
(measured by m^jm p etc.) than those with staggered quarks. This leads to the difference 
in the conclusions about the location of the physical point in Fig. 13(a). On the other 
hand, both of these studies discuss sizable deviation of several physical observables from the 
experimental values, meaning that the deviation from the continuum limit is large at N t = 4 
where these simulations are done. We should certainly carry out a calculation at larger N t 
or with an improved action in order to draw a definite conclusion about the nature of the 
QCD transition in the real world. 
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§6. Summary 

With the present power of computers, we can perform reliable extrapolations to the con- 
tinuum limit for several physical quantities in the quenched approximation of QCD. Precise 
values of the transition temperature, latent heat, interface tension, etc. obtained from de- 
tailed finite size scaling analyses are discussed in the literature. For full QCD simulations, 
however, several hundred times more computer time is required compared to quenched sim- 
ulations. This corresponds to about 5-10 years difference in the development of computer 
speed. While a few projects to construct a computer with such a speed have been proposed, 
we are not simply waiting for new computers. Many theoretical developments, especially the 
progress in improving lattice actions, open us the possibility to begin realistic simulations 
of full QCD on present supercomputers. Applications of improved actions to full QCD at 
finite temperatures have just begun. 
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